#pragma once

#include <cassert>
#include <cmath>
#include <cstdlib>
#include <ctime>
#include <initializer_list>
#include <iomanip>
#include <iostream>
#include <utility>

using Real = double;

#define SIGN(a) ((a >= 0) ? 1 : -1)
#define MAX(a, b) ((a) > (b) ? (a) : (b))
#define MIN(a, b) ((a) > (b) ? (b) : (a))

#define ACCURACY 1e-12 // set the ACCURACY
#define ITERATION 100  // set iterate times

class Complex
{
protected:
    double real, image;

public:
    Complex(double first = 0, double second = 0, bool ifPolar = false)
        : real((ifPolar ? first * cos(second) : first)), image((ifPolar ? first * sin(second) : second)) {}
    Complex(const Complex &complex) : real(complex.real), image(complex.image) {}

// 加减法
#define COMPLEX_BINARY_OP(OpName, Op)                             \
    template <class T1, class T2>                                 \
    friend Complex OpName(const T1 &lhs, const T2 &rhs)           \
    {                                                             \
        Complex z1(lhs), z2(rhs);                                 \
        return Complex(z1.real Op z2.real, z1.image Op z2.image); \
    }

    COMPLEX_BINARY_OP(operator+, +)
    COMPLEX_BINARY_OP(operator-, -)
#undef COMPLEX_BINARY_OP

// 比较操作
#define COMPLEX_COMPARE_OP(OpName, Op)               \
    template <class T1, class T2>                    \
    friend bool OpName(const T1 &lhs, const T2 &rhs) \
    {                                                \
        Complex z1(lhs), z2(rhs);                    \
        return (z1.real Op z2.real ? true : false);  \
    }

    COMPLEX_COMPARE_OP(operator>, >)
    COMPLEX_COMPARE_OP(operator<, <)
    COMPLEX_COMPARE_OP(operator>=, >=)
    COMPLEX_COMPARE_OP(operator<=, <=)
#undef COMPLEX_COMPARE_OP

    template <class T1, class T2>
    friend bool operator==(const T1 &lhs, const T2 &rhs)
    {
        Complex z1(lhs), z2(rhs);
        return ((z1.real == z2.real) && (z1.image == z2.image) ? true : false);
    }

    template <class T1, class T2>
    friend Complex operator*(const T1 &lhs, const T2 &rhs)
    {
        Complex z1(lhs), z2(rhs);
        return Complex(z1.real * z2.real - z1.image * z2.image,
                       z1.real * z2.image + z1.image * z2.real);
    }

    template <class T1, class T2>
    friend Complex operator/(const T1 &lhs, const T2 &rhs)
    {
        Complex z1(lhs), z2(rhs);
        double module = z2.real * z2.real + z2.image * z2.image;
        return Complex((z1.real * z2.real + z1.image * z2.image) / module,
                       (z1.image * z2.real - z1.real * z2.image) / module);
    }

    template <class T>
    void operator=(const T &rhs)
    {
        Complex r(rhs);
        real = r.real;
        image = r.image;
    }

    double getReal() const { return real; }
    double getImage() const { return image; }

    bool ifReal() const { return (image == 0 ? true : false); }
    bool ifImage() const { return (real == 0 ? true : false); }

    friend std::ostream &operator<<(std::ostream &out, const Complex &complex)
    {
        out << complex.real;
        if (!complex.ifReal())
        {
            out << "," << complex.image << "i";
        }
        return out;
    }

    static inline double module(const Complex &z)
    {
        return sqrt(z.getReal() * z.getReal() + z.getImage() * z.getImage());
    }

    static inline Complex normalize(const Complex &z)
    {
        return Complex(z / module(z));
    }

    static inline double arg(const Complex &z)
    {
        return atan2(z.getImage(), z.getReal());
    }

    static inline Complex conjugate(const Complex &z)
    {
        return Complex(z.getReal(), -1 * z.getImage());
    }

    static inline Complex power(const Complex &z, Real p)
    {
        return Complex(pow(module(z), p), arg(z) * p, true);
    }

    static inline auto sqrt_C(const Real &t)
    {
        return (t < 0 ? Complex(0, sqrt(-1 * t)) : sqrt(t));
    }
    static inline Complex sqrt_C(const Complex &z)
    {
        return Complex(sqrt(module(z)), arg(z) / 2, true);
    }

    static inline Complex rotate(const Complex &z, Real theta)
    {
        return Complex(module(z), arg(z) + theta, true);
    }
};

// =====================================================================================================

template <typename T>
class Node
{
public:
    T value;
    Node *next;

    Node(const T &value = T(), Node *next = nullptr) : next(next), value(value) {}

    ~Node()
    {
        if (next != nullptr)
        {
            delete next;
            next = nullptr;
        }
    }
};

template <typename T>
class Array
{
protected:
    Node<T> *head;
    int Size;

    Node<T> *previous(Node<T> *n) { return getNode(indexOf(n->value) - 1); }

    Node<T> *getNode(int i)
    {
        Node<T> *p = head;
        for (int j = 0; j < Size; j++)
        {
            if (j == i)
            {
                return p;
            }
            p = p->next;
        }
        return nullptr;
    }

    Node<T> *nodeOf(T value)
    {
        Node<T> *p = head;
        for (int i = 0; i < Size; i++)
        {
            if (p->value == value)
            {
                return p;
            }
            p = p->next;
        }
        return p;
    }

public:
    Array() : head(nullptr), Size(0) {}
    Array(const int Size, const T &t = T()) : Size(Size)
    {
        head = new Node<T>(t);
        Node<T> *p = head;
        for (int i = 1; i < Size; i++)
        {
            p->next = new Node<T>(t);
            p = p->next;
        }
    }
    Array(std::initializer_list<T> l) : Size(l.size())
    {
        auto j = l.begin();
        head = new Node<T>(*j++);
        Node<T> *p = head;
        for (int i = 1; i < Size; i++)
        {
            p->next = new Node<T>(*j++);
            p = p->next;
        }
    }
    Array(const Array &arr) : Size(arr.size())
    {
        head = (arr.empty() ? nullptr : new Node<T>(arr[0]));
        Node<T> *p = head;
        for (int i = 1; i < Size; i++)
        {
            p->next = new Node<T>(arr[i]);
            p = p->next;
        }
    }

    T &operator[](int i)
    {
        Node<T> *p = head;
        while (i > 0)
        {
            if (p->next == nullptr)
            {
                p->next = new Node<T>;
                Size++;
            }
            p = p->next;
            i--;
        }
        return p->value;
    }
    T operator[](int i) const
    {
        Node<T> *p = head;
        while (i > 0)
        {
            p = p->next;
            i--;
        }
        return p->value;
    }

    void operator=(const Array &arr)
    {
        Array tmp(arr);
        head = tmp.head;
    }

    bool operator==(const Array &arr)
    {
        Node<T> *p = head;
        if (Size != arr.size())
        {
            return false;
        }
        for (int i = 0; i < Size; i++)
        {
            if ((*this)[i] != arr[i])
            {
                return false;
            }
        }
        return true;
    }

    friend std::ostream &operator<<(std::ostream &out, const Array<T> &arr)
    {
        out << "[";
        for (int i = 0; i < arr.size(); i++)
        {
            out << arr[i] << " ";
        }
        out << "\b]";
        return out;
    }

    void insert(int i, T value)
    {
        Node<T> *q = new Node<T>(value);
        if (i == 0)
        {
            q->next = head;
            head = q;
        }
        else
        {
            Node<T> *p = getNode(i - 1);
            q->next = p->next;
            p->next = q;
        }
        Size++;
    }

    int indexOf(T value)
    {
        for (int i = 0; i < Size; i++)
        {
            if ((*this)[i] == value)
            {
                return i;
            }
        }
        return -1;
    }

    T pop()
    {
        if (Size == 0)
        {
            return 0;
        }
        if (Size == 1)
        {
            delete head;
            head = nullptr;
        }
        Node<T> *p = getNode(Size - 2);
        T value = p->next->value;
        delete p->next;
        p->next = nullptr;
        Size--;
        return value;
    }

    void append(T value)
    {
        Node<T> *q = new Node<T>(value);
        if (Size == 0)
        {
            head = q;
        }
        else
        {
            Node<T> *p = getNode(Size - 1);
            p->next = q;
        }
        Size++;
    }

    bool kill(T value)
    {
        if (Size == 0)
        {
            return false;
        }
        Node<T> *p = nodeOf(value);
        if (p == nullptr)
        {
            return false;
        }
        Node<T> *q = previous(p);
        (q == nullptr ? head = head->next : q->next = p->next);
        p->next = nullptr;
        delete p;
        p = nullptr;
        Size--;
        return true;
    }

    bool remove(int i)
    {
        if (i + 1 > Size)
        {
            return false;
        }
        Node<T> *tmp, *p;
        if (i == 0)
        {
            tmp = head;
            head = head->next;
        }
        else
        {
            p = getNode(i - 1);
            tmp = p->next;
            p->next = tmp->next;
        }
        tmp->next = nullptr;
        delete tmp;
        tmp = nullptr;
        Size--;
        return true;
    }

    int size() const { return Size; }
    bool empty() const { return (head == nullptr) ? true : false; }

    void clean()
    {
        if (head->next != nullptr)
        {
            delete head->next;
            head->value = 0;
        }
    }

    ~Array()
    {
        delete head;
        head = nullptr;
    }
};

// ================================================================================================

/**
 * @file   Matrix.h
 * @author xingyifan <xingyifan@xingyifan-VirtualBox>
 * @date   Wed Jul 14 10:34:03 2021
 *
 * @brief  vector
 *
 *
 */

template <class T>
class Vector
{
protected:
    const int Dim;
    T *value;

public:
    // 初始化构造函数
    Vector(const int Dim, const T &t = T()) : Dim(Dim), value(new T[Dim])
    {
        for (int i = 0; i < Dim; value[i++] = t)
            ;
    }

    // 数组初始化
    Vector(const int Dim, const T arr[]) : Dim(Dim), value(new T[Dim])
    {
        for (int i = 0; i < Dim; value[i] = arr[i], i++)
            ;
    }

    // 可以通过列表初始化
    Vector(std::initializer_list<T> l) : Dim(l.size()), value(new T[Dim])
    {
        auto j = l.begin();
        for (int d = 0; d < Dim; ++d)
            value[d] = *j++;
    }

    // 一定要注意还有普通拷贝构造函数！！！！
    Vector(const Vector &rhs) : Dim(rhs.size()), value(new T[Dim])
    {
        for (int d = 0; d < Dim; d++)
            value[d] = rhs[d];
    }

    // 只能显式地调用此构造函数
    template <class T2>
    explicit Vector(const Vector<T2> &rhs) : Dim(rhs.size()), value(new T[Dim])
    {
        for (int d = 0; d < Dim; d++)
            value[d] = static_cast<T>(rhs[d]);
    }

    // 构造第 D 位为 1 的单位向量
    static Vector<T> unit(int Dim, int D = 0)
    {
        Vector<T> r(Dim);
        r[D] = static_cast<T>(1);
        return r;
    }

    T &operator[](int i) { return value[i]; }
    T operator[](int i) const { return value[i]; }

// 加一个向量
#define VECTOR_BINARY_OP(OpName, Op)             \
    template <class T2>                          \
    auto OpName(const Vector<T2> &rhs) const     \
    {                                            \
        using Tx = decltype(value[0] Op rhs[0]); \
        Vector<Tx> res(Dim);                     \
        for (int i = 0; i < Dim; i++)            \
            res[i] = value[i] Op rhs[i];         \
        return res;                              \
    }

    VECTOR_BINARY_OP(operator+, +)
    VECTOR_BINARY_OP(operator-, -)
    VECTOR_BINARY_OP(operator*, *)
    VECTOR_BINARY_OP(operator/, /)
#undef VECTOR_BINARY_OP

// 加一个常数
#define RIGHT_BROADCAST(OpName, Op)           \
    template <class T2>                       \
    auto OpName(const T2 &rhs) const          \
    {                                         \
        using Tx = decltype(value[0] Op rhs); \
        Vector<Tx> res(Dim);                  \
        for (int d = 0; d < Dim; ++d)         \
            res[d] = value[d] Op rhs;         \
        return res;                           \
    }

    RIGHT_BROADCAST(operator+, +)
    RIGHT_BROADCAST(operator-, -)
    RIGHT_BROADCAST(operator*, *)
    RIGHT_BROADCAST(operator/, /)
#undef RIGHT_BROADCAST

    friend std::ostream &operator<<(std::ostream &out, const Vector<T> &vec)
    {
        for (int i = 0; i < vec.size(); i++)
        {
            out << vec[i] << " ";
        }
        out << std::endl;
        return out;
    }

    void operator=(const Vector &rhs)
    {
        for (int d = 0; d < Dim; d++)
            value[d] = rhs[d];
    }

    // 获取一定范围的向量切片
    Vector<T> range(int left, int right)
    {
        Vector<T> res(right - left);
        for (int i = left; i < right; i++)
        {
            res[i - left] = value[i];
        }
        return res;
    }

    // 检查是否在向量中
    bool find(T t) const
    {
        for (int i = 0; i < Dim; i++)
        {
            if (t == value[i])
            {
                return true;
            }
        }
        return false;
    }

    int size() const { return Dim; }

    ~Vector() { delete[] value; }
};

// ===================================================================================================

// 置换向量（每一位存放所在的位，例如 p[0] = 5，表示 5 元素在 0 位置）
class Perm : public Vector<int>
{
public:
    // 初始化构造函数
    Perm(const int Dim) : Vector<int>(Dim)
    {
        for (int i = -1; i < Dim - 1; this->value[++i] = i)
            ;
    }

    // 拷贝构造
    Perm(const Perm &perm) : Vector<int>(perm.size())
    {
        for (int i = -1; i < Dim - 1; this->value[++i] = perm[i])
            ;
    }

    // 对换
    void exchange(int i, int j)
    {
        int tmp = this->value[i];
        this->value[i] = this->value[j];
        this->value[j] = tmp;
    }

    // 置换的逆
    static Perm inverse(const Perm &p)
    {
        int Dim = p.size();
        Perm q(Dim);
        for (int i = 0; i < Dim; i++)
        {
            q[p[i]] = i;
        }
        return q;
    }

// 置换作用
#define PERM_BINARY_OP(Type, OpName, Op)    \
    auto OpName(const Type &rhs) const      \
    {                                       \
        Type p_rhs(this->Dim);              \
        for (int i = 0; i < this->Dim; i++) \
        {                                   \
            p_rhs[i] = rhs[this->value[i]]; \
        }                                   \
        return p_rhs;                       \
    }

    template <class T>
    PERM_BINARY_OP(Vector<T>, operator*, *)
    PERM_BINARY_OP(Perm, operator*, *)
#undef PERM_BINARY_OP
};

// 如果要在头文件中写函数，最好写成内联函数，避免重复定义
struct Vec
{
    template <class T>
    static bool approx(const Vector<T> &lhs, const Vector<T> &rhs)
    {
        if (rhs.size() != lhs.size())
        {
            return false;
        }
        for (int i = 0; i < lhs.size(); i++)
        {
            if ((lhs[i] - rhs[i]) > ACCURACY)
            {
                return false;
            }
        }
        return true;
    }

    static Vector<Real> rand_vec(int Dim)
    {
        Vector<Real> r(Dim);
        srand((int)time(0));
        for (int i = 0; i < Dim; i++)
        {
            r[i] = rand() / Real(RAND_MAX);
        }
        return r;
    }

    template <class T, int Dim>
    static auto cross(const Vector<T> &va, const Vector<T> &vb)
    {
        if (Dim == 2)
        {
            return va[0] * vb[1] - va[1] * vb[0];
        }
        else if (Dim == 3)
        {
            return Vector<T>{
                va[1] * vb[2] - va[2] * vb[1],
                va[2] * vb[0] - va[0] * vb[2],
                va[0] * vb[1] - va[1] * vb[0]};
        }
        return 0;
    }

    // 幂
    template <class T>
    static Vector<T> power(const Vector<T> &rhs, Real p)
    {
        int Dim = rhs.size();
        Vector<T> res(Dim);
        for (int i = 0; i < Dim; res[i] = pow(static_cast<Real>(rhs[i]), p), i++)
            ;
        return res;
    }

    // 去掉右边的零进行压缩
    template <class T>
    static Vector<T> nonzero(const Vector<T> &rhs)
    {
        int Dim = rhs.size();
        for (int i = Dim - 1; i >= 0; i--)
        {
            if (rhs[i] != 0)
            {
                return rhs.range(0, i + 1);
            }
        }
        return rhs;
    }

    template <class T>
    static Vector<T> reverse(const Vector<T> &rhs)
    {
        int Dim = rhs.size();
        Vector<T> res(Dim);
        for (int i = 0; i < Dim; res[i] = rhs[Dim - i - 1], i++)
            ;
        return res;
    }

    template <class T>
    static T sum(const Vector<T> &rhs)
    {
        T a = 0;
        for (int i = 0; i < rhs.size(); i++)
        {
            a += rhs[i];
        }
        return a;
    }

    template <class T>
    static T dot(const Vector<T> &vec1, const Vector<T> &vec2, Real weight = 1)
    {
        int Dim = vec1.size();
        T sum = 0;
        for (int i = 0; i < Dim; i++)
        {
            sum += vec1[i] * vec2[i] * weight;
        }
        return sum;
    }

    template <class T>
    static Vector<int> sgn(const Vector<T> &rhs)
    {
        int Dim = rhs.size();
        Vector<int> res(Dim);
        for (int i = 0; i < Dim; i++)
        {
            res[i] = (rhs[i] > 0 ? 1 : (rhs[i] < 0 ? -1 : 0));
        }
        return res;
    }

    template <class T>
    static Vector<int> abs(const Vector<T> &rhs)
    {
        int Dim = rhs.size();
        Vector<int> res(Dim);
        for (int d = 0; d < Dim; res[d] = fabs(rhs[d]), ++d)
            ;
        return res;
    }

    template <class T>
    static T norm(const Vector<T> &lhs, int nt = 2)
    {
        int Dim = lhs.size();
        if (nt == 2)
        {
            return sqrt(dot(lhs, lhs));
        }
        else if (nt == 1)
        {
            T a = 0;
            for (int i = 0; i < Dim; a += fabs(lhs[i]), ++i)
                ;
            return a;
        }
        else if (nt == 0)
        {
            T a = fabs(lhs[0]);
            for (int i = 1; i < Dim; i++)
            {
                a = MAX(fabs(lhs[i]), a);
            }
            return a;
        }
        return 0; // never reached
    }

    template <class T>
    static Vector<T> normalize(const Vector<T> &rhs, int nt = 2)
    {
        T normalizer = norm(rhs, nt);
        return Vector<T>(rhs / normalizer);
    }

    template <class T, class Cmp = std::less<T>>
    static T min_of(const Vector<T> &v, Cmp cmp = Cmp())
    {
        int Dim = v.size();
        auto a = v[0];
        for (int d = 1; d < Dim; d++)
            if (cmp(v[d], a))
                a = v[d];
        return a;
    }

    template <class T>
    static T max_of(const Vector<T> &v)
    {
        return min_of(v, std::greater<T>());
    }

    template <class T>
    static int max_index_of(const Vector<T> &v)
    {
        int Dim = v.size();
        auto a = v[0];
        int index = 0;
        for (int d = 1; d < Dim; d++)
        {
            index = (fabs(v[d]) > fabs(v[index]) ? d : index);
        }
        return index;
    }

    // element-wise min
    template <class T>
    static Vector<T> min(const Vector<T> &lhs, const Vector<T> &rhs)
    {
        int Dim = lhs.size();
        Vector<T> res(Dim);
        for (int d = 0; d < Dim; d++)
        {
            res[d] = MIN(lhs[d], rhs[d]);
        }
        return res;
    }

    // element-wise max
    template <class T>
    static Vector<T> max(const Vector<T> &lhs, const Vector<T> &rhs)
    {
        int Dim = lhs.size();
        Vector<T> res(Dim);
        for (int d = 0; d < Dim; d++)
        {
            res[d] = MAX(lhs[d], rhs[d]);
        }
        return res;
    }
};

// =============================================================================================

template <class T>
class Matrix : public Vector<T>
{
protected:
    const int Row, Col;

public:
    Matrix(const int Row, const int Col, const T &t = T()) : Row(Row), Col(Col), Vector<T>(Row * Col, t) {} // 初始化构造函数
    Matrix(const int Row, const int Col, std::initializer_list<T> l) : Row(Row), Col(Col), Vector<T>(l) {}  // 通过列表初始化
    Matrix(const int Row, const int Col, const T arr[]) : Row(Row), Col(Col), Vector<T>(Row * Col, arr) {}  // 通过列表初始化

    // 一定要注意还有普通拷贝构造
    Matrix(const Matrix &rhs) : Row(rhs.row()), Col(rhs.col()), Vector<T>(rhs.size())
    {
        for (int d = 0; d < this->Dim; d++)
            this->value[d] = rhs[d];
    }

    // 只能显式地调用此构造函数
    template <class T2>
    explicit Matrix(const Matrix<T2> &rhs) : Row(rhs.row()), Col(rhs.col()), Vector<T>(rhs.size())
    {
        for (int d = 0; d < this->Dim; d++)
            this->value[d] = static_cast<T>(rhs[d]);
    }

    // 返回单位阵
    static Matrix<T> unit(int row)
    {
        Matrix<T> I(row, row);
        for (int d = 0; d < row; d++)
            I(d, d) = 1;
        return I;
    }

    // 注意在模板类中无法直接访问基类成员，需要通过 this 指针调用
    T &operator()(int i, int j) { return this->value[i * Row + j]; }
    T operator()(int i, int j) const { return this->value[i * Row + j]; }

    // 返回第 j 列
    Vector<T> operator()(int j) const
    {
        Vector<T> v(Row);
        for (int i = 0; i < Row; i++)
        {
            v[i] = this->value[i * Row + j];
        }
        return v;
    }

    // 在 r,c 位置插入列向量 v
    template <class T2>
    void insert(int r, int c, Vector<T2> &v)
    {
        int maxRow = MIN(v.size(), Row - r);
        for (int i = r; i < maxRow + r; i++)
        {
            this->value[r * Row + c] = v[i - r];
        }
    }

    // 获取矩阵的子阵
    Matrix<T> subMatirx(Vector<int> &horizon, Vector<int> &vertical)
    {
        int rowNum = horizon.size(), colNum = vertical.size();
        Matrix<T> sub(rowNum, colNum);
        for (int i = 0; i < rowNum; i++)
        {
            for (int j = 0; j < colNum; j++)
            {
                sub(i, j) = this->value[horizon[i] * Row + vertical[j]];
            }
        }
        return sub;
    }

    int row() const { return Row; }
    int col() const { return Col; }
    int order() const { return MIN(Row, Col); }
    bool square() const { return (Row == Col ? true : false); }

    friend std::ostream &operator<<(std::ostream &out, const Matrix<T> &mat)
    {
        int Row = mat.row(), Col = mat.col();
        for (int i = 0; i < Row; i++)
        {
            for (int j = 0; j < Col; j++)
            {
                out << std::setw(12) << mat(i, j) << " ";
            }
            out << std::endl;
        }
        return out;
    }

// 加一个矩阵/在对角线上加一个向量
#define MATRIX_BINARY_OP(OpName, Op)                       \
    template <class T2>                                    \
    auto OpName(const Matrix<T2> &rhs) const               \
    {                                                      \
        using Tx = decltype(this->value[0] Op rhs[0]);     \
        Matrix<Tx> res(Row, Col);                          \
        for (int i = 0; i < this->Dim; i++)                \
            res[i] = this->value[i] Op rhs[i];             \
        return res;                                        \
    }                                                      \
    template <class T2>                                    \
    auto OpName(const Vector<T2> &rhs) const               \
    {                                                      \
        using Tx = decltype(this->value[0] Op rhs[0]);     \
        Matrix<Tx> res(Row, Col, this->value);             \
        for (int i = 0; i < rhs.size(); i++)               \
            res[i * Row + i] = res[i * Row + i] Op rhs[i]; \
        return res;                                        \
    }

    MATRIX_BINARY_OP(operator+, +)
    MATRIX_BINARY_OP(operator-, -)
#undef MATRIX_BINARY_OP

// 矩阵赋值（定义矩阵以后，可以直接用向量赋值）
#define MATRIX_EVALUATE(Type)               \
    void operator=(const Type &rhs)         \
    {                                       \
        for (int i = 0; i < this->Dim; i++) \
        {                                   \
            this->value[i] = rhs[i];        \
        }                                   \
    }

    MATRIX_EVALUATE(Vector<T>)
    MATRIX_EVALUATE(Matrix)
#undef MATRIX_EVALUATE

    // 乘除法
#define MATRIX_BINARY_OP(OpName, Op)        \
    auto OpName(const T &rhs) const         \
    {                                       \
        Matrix res(Row, Col);               \
        for (int i = 0; i < this->Dim; i++) \
            res[i] = this->value[i] Op rhs; \
        return res;                         \
    }

    MATRIX_BINARY_OP(operator*, *)
    MATRIX_BINARY_OP(operator/, /)
#undef MATRIX_BINARY_OP

    // 矩阵乘法
    Vector<T> operator*(const Vector<T> &_v)
    {
        Vector<T> v(Row);
        T sum;
        for (int i = 0; i < Row; i++)
        {
            sum = 0;
            for (int j = 0; j < Col; j++)
            {
                sum += operator()(i, j) * _v[j];
            }
            v[i] = sum;
        }
        return v;
    }
    Matrix operator*(const Matrix &_m)
    {
        int _column = _m.col();
        Matrix m(Row, _column);
        T sum;

        for (int i = 0; i < Row; i++)
        {
            for (int j = 0; j < _column; j++)
            {
                sum = 0;
                for (int k = 0; k < Col; k++)
                {
                    sum += operator()(i, k) * _m(k, j);
                }
                m(i, j) = sum;
            }
        }
        return m;
    }
};

// ====================================================================================================

struct Mat
{
    template <class T>
    static bool approx(const Matrix<T> &lhs, const Matrix<T> &rhs)
    {
        if (rhs.row() != lhs.row() || rhs.col() != lhs.col())
        {
            return false;
        }
        for (int i = 0; i < lhs.row(); i++)
        {
            for (int j = 0; j < lhs.col(); j++)
            {
                if ((lhs(i, j) - rhs(i, j)) > ACCURACY)
                {
                    return false;
                }
            }
        }
        return true;
    }

    // 置换作用于矩阵
    template <class T>
    static Matrix<T> exchange(const Perm &perm, const Matrix<T> &A)
    {
        int row = A.row(), col = A.col();
        Matrix<T> pA(row, col);
        for (int i = 0; i < row; i++)
        {
            for (int j = 0; j < col; j++)
            {
                pA(i, j) = A(perm[i], j);
            }
        }
        return pA;
    }

    static Matrix<Real> Hilbert(int n)
    {
        Matrix<Real> A(n, n);
        for (int i = 0; i < n; i++)
        {
            for (int j = 0; j < n; j++)
            {
                A(i, j) = 1.0 / (i + j + 1);
            }
        }
        return A;
    }

    // 多项式的拉普拉斯阵
    static Matrix<Real> Laplace(int n, Real a[])
    {
        Matrix<Real> A(n, n);
        for (int i = 0; i < n; i++)
        {
            A(0, i) = a[i];
            if (i < n - 1)
            {
                A(i + 1, i) = -1;
            }
        }
        return A;
    }

    static Matrix<Real> rand_mat(int row, int column)
    {
        Matrix<Real> m(row, column);
        srand((int)time(0));
        for (int i = 0; i < row; i++)
        {
            for (int j = 0; j < column; j++)
            {
                m(i, j) = rand() / Real(RAND_MAX);
            }
        }
        return m;
    }

    // ====================================================================================

    template <class T>
    static T norm(const Matrix<T> &lhs, int nt = 2)
    {
        int Row = lhs.row(), Col = lhs.col();
        if (nt == 2)
        {
            return sqrt(Vec::dot(lhs, lhs));
        }
        else if (nt == 1)
        {
            T max = 0;
            for (int i = 0; i < Col; i++)
            {
                T a = 0;
                for (int j = 0; j < Row; a += fabs(lhs(j++, i)))
                    ;
                max = (a > max ? a : max);
            }
            return max;
        }
        else if (nt == 0)
        {
            T max = 0;
            for (int i = 0; i < Row; i++)
            {
                T a = 0;
                for (int j = 0; j < Col; a += fabs(lhs(i, j++)))
                    ;
                max = (a > max ? a : max);
            }
            return max;
        }
        return 0; // never reached
    }

    // e 范数 / f 范数
    template <class T>
    static T norm(const Matrix<T> &lhs, char c)
    {
        int Row = lhs.row(), Col = lhs.col();
        if (c == 'e')
        {
            T sum = 0;
            for (int i = 0; i < Row; i++)
            {
                for (int j = 0; j < Col; j++)
                {
                    if (i != j)
                    {
                        sum += lhs(i, j) * lhs(i, j);
                    }
                }
            }
            return sqrt(sum);
        }
        else if (c == 'f')
        {
            T sum = 0;
            for (int i = 0; i < Row; i++)
            {
                for (int j = 0; j < Col; j++)
                {
                    sum += lhs(i, j) * lhs(i, j);
                }
            }
            return sqrt(sum);
        }
        return 0;
    }

    // get upper or lower
    template <class T>
    static Matrix<T> upper(Matrix<T> &A, bool hasDialog = true)
    {
        int Col = A.col(), Row = A.col();
        Matrix<T> m(Row, Col);
        int limit = MIN(Row, Col);
        for (int i = 0; i < limit; i++)
        {
            for (int j = (hasDialog ? i : i + 1); j < limit; j++)
            {
                m(i, j) = A(i, j);
            }
        }
        return m;
    }

    template <class T>
    static Matrix<T> lower(Matrix<T> &A, bool hasDialog = true)
    {
        int Col = A.col(), Row = A.col();
        Matrix<T> m(Row, Col);
        int limit = MIN(Row, Col);
        for (int i = 0; i < limit; i++)
        {
            for (int j = 0; j < (hasDialog ? i + 1 : i); j++)
            {
                m(i, j) = A(i, j);
            }
        }
        return m;
    }

    template <class T>
    static Vector<T> diag(Matrix<T> &A)
    {
        int Col = A.col(), Row = A.col();
        int limit = MIN(Row, Col);
        Vector<T> v(limit);
        for (int i = 0; i < limit; i++)
        {
            v[i] = A(i, i);
        }
        return v;
    }

    template <class T>
    static Matrix<T> transpose(Matrix<T> &A)
    {
        int Col = A.col(), Row = A.col();
        Matrix<T> m(Col, Row);
        for (int i = 0; i < Col; i++)
        {
            for (int j = 0; j < Row; j++)
            {
                m(i, j) = A(j, i);
            }
        }
        return m;
    }

    template <class T>
    static Matrix<T> power(Matrix<T> &A, int p)
    {
        Matrix<T> P = Matrix<T>::unit(A.row());
        for (int i = 0; i < p; i++)
        {
            P = P * A;
        }
        return P;
    }

    template <class T>
    static T det(Matrix<T> &A)
    {
        int sgn = -1, n = A.row();
        if (n == 1)
        {
            return A[0];
        }
        T sum = 0;
        Vector<int> rowIndex(n - 1), colIndex(n - 1);
        // 子矩阵的行列向量标记
        for (int i = 1; i < n; i++)
        {
            rowIndex[i - 1] = i;
            colIndex[i - 1] = i;
        }
        Matrix<T> B(n - 1, n - 1);
        for (int i = 0; i < n; i++)
        {
            sgn *= -1;
            B = A.subMatirx(rowIndex, colIndex);
            sum += sgn * A[i] * det(B);
            colIndex[i] = i;
        }
        return sum;
    }

    template <class T>
    Complex eigen2(const Matrix<T> &A, int row = 0, int col = 0)
    {
        T b = -1 * (A(row, col) + A(row + 1, col + 1));
        T c = A(row, col) * A(row + 1, col + 1) - A(row, col + 1) * A(row + 1, col);
        Complex z = sqrt_C(b * b - 4 * c);
        z = (z - b) / 2;
        return z;
    }
};

// =========================================================================================

/**
 * @file   Matrix.h
 * @author xingyifan <xingyifan@xingyifan-VirtualBox>
 * @date   Wed Jul 14 10:33:43 2021
 *
 * @brief  sparse matrix
 *
 *
 */

class SpMatrix
{
private:
    const int m, n;
    int nnz;
    Array<int> IA, JA;
    Array<double> AA; //对角元放在每行首位

public:
    SpMatrix(int row, int column) : m(row), n(column)
    {
        nnz = std::min(m, n);
        for (int i = 0; i < m; i++)
        {
            IA[i] = i;
        }
        for (int i = 0; i < nnz; i++)
        {
            JA[i] = i;
            AA[i] = 0;
        }
        JA[nnz] = __INT_MAX__;
        IA[m] = nnz;
        AA[nnz] = __FLT_MAX__;
    }
    SpMatrix(int _row, int _column, int _nnz, int *_IA, int *_JA, double *_AA) : m(_row), n(_column)
    {
        nnz = _nnz;
        for (int i = 0; i < m; i++)
        {
            IA[i] = _IA[i];
        }
        IA[m] = nnz;
        for (int i = 0; i < n; i++)
        {
            JA[i] = _JA[i];
        }
        JA[nnz] = __INT_MAX__;
        for (int i = 0; i < nnz; i++)
        {
            AA[i] = _AA[i];
        }
        AA[nnz] = __FLT_MAX__;
    }
    SpMatrix(const SpMatrix &spm) : m(spm.m), n(spm.n)
    {
        nnz = spm.nnz;
        IA = spm.IA;
        JA = spm.JA;
        AA = spm.AA;
    }

// 加减
#define SPMATRIX_BINARY_OP(OpName, Op)                            \
    auto OpName(const SpMatrix &_spm) const                       \
    {                                                             \
        SpMatrix spm = *this;                                     \
        for (int i = 0; i < m; i++)                               \
        {                                                         \
            int range = _spm.IA[i + 1] - _spm.IA[i];              \
            for (int j = 0; j < range; j++)                       \
            {                                                     \
                int col = _spm.IA[i] + j;                         \
                spm.insert(i, _spm.JA[col], Op 1 * _spm.AA[col]); \
            }                                                     \
        }                                                         \
        return spm;                                               \
    }

    SPMATRIX_BINARY_OP(operator+, +)
    SPMATRIX_BINARY_OP(operator-, -)
#undef SPMATRIX_BINARY_OP

    SpMatrix operator*(double k)
    {
        SpMatrix spm = *this;
        for (int i = 0; i < nnz; i++)
        {
            spm.AA[i] *= k;
        }
        return spm;
    }
    //稀疏矩阵的乘法
    SpMatrix operator*(const SpMatrix &_spm)
    {
        const int _m = m;
        const int _n = _spm.n;
        SpMatrix spm(_m, _n);
        for (int i = 0; i < m; i++)
        {
            double tmp_val[_n];
            //计算i行的结果
            for (int j = IA[i]; j < IA[i + 1]; j++)
            {
                //这是一个很巧妙的计算方式，只需要每次计算两个元素的乘积，记录在对应的列的位置上
                //然后不断累加即可，相当于将n次乘法和加法分开计算
                for (int k = _spm.IA[_spm.JA[j]]; k < _spm.IA[_spm.JA[j] + 1]; k++)
                {
                    tmp_val[_spm.JA[k]] += AA[j] * _spm.AA[k];
                }
            }
            //去掉过小的元素
            for (int j = 0; j < _n; j++)
            {
                if (fabs(tmp_val[j]) > __FLT_MIN__)
                {
                    //插入结果
                    spm.insert(i, j, tmp_val[j]);
                }
            }
        }
        return spm;
    }

    friend std::ostream &operator<<(std::ostream &out, const SpMatrix &_spm)
    {
        for (int i = 0; i < _spm.m; i++)
        {
            // IA[m] = nnz，边界值
            for (int j = _spm.IA[i]; j < _spm.IA[i + 1]; j++)
            {
                out << i << ", " << _spm.JA[j] << ", " << _spm.AA[j] << std::endl;
            }
        }
        return out;
    }

    double get(int i, int j) const
    {
        //只能获取存在的位置元素
        for (int k = IA[i]; k < IA[i + 1]; k++)
        {
            if (JA[k] == j)
            {
                return AA[k];
            }
        }
        return 0;
    }

    double &set(int i, int j)
    {
        //只能设置存在的位置
        for (int k = IA[i]; k < IA[i + 1]; k++)
        {
            if (JA[k] == j)
            {
                return AA[k];
            }
        }
        return AA[nnz];
    }

    void insert(int i, int j, double value)
    {
        //如果插入到存在的元素位置，则直接加上插入值
        if (i == j)
        {
            int k = IA[i];
            AA[k] += value;
            return;
        }
        int ins = IA[i] + 1;
        for (int k = ins; k < IA[i + 1]; k++)
        {
            if (JA[k] == j)
            {
                AA[k] += value;
                return;
            }
            else if (JA[k] > j)
            {
                //利用动态数组直接插入
                JA.insert(k, j);
                AA.insert(k, value);
                //改变IA和nnz
                for (int p = i + 1; p < m + 1; p++)
                {
                    IA[p]++;
                }
                nnz++;
                return;
            }
        }
        //当只有对角元的时候插入到下一行开头位置
        JA.insert(ins, j);
        AA.insert(ins, value);
        //改变IA和nnz
        for (int p = i + 1; p < m + 1; p++)
        {
            IA[p]++;
        }
        nnz++;
    }

    SpMatrix transpose()
    {
        SpMatrix res(n, m);
        for (int i = 0; i < m; i++)
        {
            for (int j = frontOfRow(i); j < frontOfRow(i + 1); j++)
            {
                res.insert(posOfCol(j), i, global(j));
            }
        }
        return res;
    }

    int row() const { return m; }
    int col() const { return n; }
    int nonzero() const { return nnz; }

    double diag(int i) { return AA[IA[i]]; }      //获取对角元
    int frontOfRow(int i) const { return IA[i]; } //获取指定行的开头
    int posOfCol(int j) const { return JA[j]; }   //通过全局指标获取元素所在的列位置
    double global(int j) const { return AA[j]; }  //通过全局指标获取元素值

    //获取(i,j)元素的全局列标
    int globalCol(int i, int j) const
    {
        for (int k = IA[i]; k < IA[i + 1]; k++)
        {
            if (JA[k] == j)
            {
                return k;
            }
        }
        return -1;
    }

    // 迭代法
    static Vector<double> Jacobi(SpMatrix &A, Vector<double> &b)
    {
        int n = A.row();
        Vector<double> x(n, 1), y(n);
        double dx;
        do
        {
            dx = 0;
            for (int i = 0; i < n; i++)
            {
                double sum = 0;
                int s_row = A.frontOfRow(i);
                //只需要计算一次加法b[i] / d[i,i]
                sum += b[i] / A.global(s_row);
                for (int j = s_row + 1; j < A.frontOfRow(i + 1); j++)
                {
                    //累计D1(L+U) * x
                    int col = A.posOfCol(j);
                    sum -= A.global(j) * x[col] / A.global(s_row);
                }
                y[i] = sum;
            }
            dx = Vec::norm(y - x, 0);
            x = y;
        } while (dx >= ACCURACY);
        return x;
    }

    static Vector<double> G_S(SpMatrix &A, Vector<double> &b)
    {
        int n = A.row();
        Vector<double> x(n, 1), y(n);
        double dx;
        do
        {
            dx = 0;
            y = x;
            for (int i = 0; i < n; i++)
            {
                double sum = 0;
                int s_row = A.frontOfRow(i);
                //只需要计算一次加法b[i] / d[i,i]
                sum += b[i] / A.global(s_row);
                for (int j = s_row + 1; j < A.frontOfRow(i + 1); j++)
                {
                    //累计D1(L+U) * x
                    int col = A.posOfCol(j);
                    sum -= A.global(j) * x[col] / A.global(s_row);
                }
                // G_S迭代法每次直接更新计算出的x[i]值
                x[i] = sum;
            }
            dx = Vec::norm(y - x, 0);
        } while (dx >= ACCURACY);
        return x;
    }

    static Vector<double> SOR(SpMatrix &A, Vector<double> &b, double w)
    {
        int n = A.row();
        Vector<double> x(n, 1), y(n);
        double dx;
        do
        {
            dx = 0;
            y = x;
            for (int i = 0; i < n; i++)
            {
                double sum = 0;
                int s_row = A.frontOfRow(i);
                //只需要计算一次加法b[i] / d[i,i]
                sum += b[i] / A.global(s_row);
                for (int j = s_row + 1; j < A.frontOfRow(i + 1); j++)
                {
                    //累计D1(L+U) * x
                    int col = A.posOfCol(j);
                    sum -= A.global(j) * x[col] / A.global(s_row);
                }
                // SOR迭代法每次计算出的x[i]值，然后直接乘w修正delta x
                x[i] = sum;
            }
            x = y + (x - y) * w;
            dx = Vec::norm(y - x, 0);
        } while (dx >= ACCURACY);
        return x;
    }
};

// ====================================================================================================

struct Solver
{
    // 矩阵操作
    static Matrix<Real> inverse(Matrix<Real> &A);
    static Real Householder(Vector<Real> &x, Vector<Real> &v);
    static Vector<Real> Givens(Vector<Real> &v, int i = 0, int j = 1);

    // 矩阵分解
    static bool LU(Matrix<Real> &A);                                                    // LU 分解
    static void LU(Matrix<Real> &A, Perm &perm);                                        // 列主元 LU 分解
    static bool Cholesky(Matrix<Real> &A);                                              // Cholesky
    static bool Cholesky_D(Matrix<Real> &A);                                            // 优化的 Cholesky
    static void QR(Matrix<Real> &A, Vector<Real> &d);                                   // QR 分解
    static void rotate(Matrix<Real> &A, Matrix<Real> &Q, int p, int q);                 // Jacobi 平面旋转
    static void HAH(Matrix<Real> &A, int sR, int eR, int sC, int eC, Vector<Real> &v);  // 标准 HAH ,可选累积正交变换
    static void HAH(Matrix<Real> &A, int sR, int eR, Vector<Real> &v);                  // 对称形式的 HAH ,可选累积正交变换
    static void Hessenberg(Matrix<Real> &A, Matrix<Real> &Q, bool ifSymmetric = false); // 上 Hessenberg 化
    static void Real_QR(Matrix<Real> &A, int l, int m);                                 // 对指定区域的双重步 QR 迭代法
    static void Wilkinson_QR(Matrix<Real> &A, int l, int m);                            // 使用 Wilkinson 位移的对称 QR 迭代

    // 常规解法
    static Vector<Real> solve_upper(Matrix<Real> &A, Vector<Real> &b, bool unit = false); // 是否按照单位上三角阵求解
    static Vector<Real> solve_lower(Matrix<Real> &A, Vector<Real> &b, bool unit = false); // 是否按照单位下三角阵求解
    static Vector<Real> solve(Matrix<Real> &A, Vector<Real> &b, bool choose = false);     // 是否选主元求解
    static Vector<Real> cholesky(Matrix<Real> &A, Vector<Real> &b, bool improve = false); // 是否优化为 LDL 分解

    // 最小二乘问题
    static Vector<Real> qr(Matrix<Real> &A, Vector<Real> &b); // QR 分解法
    static Vector<Real> cg(Matrix<Real> &A, Vector<Real> &b); // 共轭梯度法

    // 迭代法
    static Vector<Real> Jacobi(Matrix<Real> &A, Vector<Real> &b);
    static Vector<Real> G_S(Matrix<Real> &A, Vector<Real> &b);
    static Vector<Real> SOR(Matrix<Real> &A, Vector<Real> &b, Real w);

    // 条件数
    static Real kappa(Matrix<Real> &A);                                               // 盲人爬山法
    static Real caculate_analysis(Matrix<Real> &A, Vector<Real> &x, Vector<Real> &b); // 线性方程的计算误差界

    // 特征值与特征向量
    static Real power(Matrix<Real> &A, Real mu = 0);                   // 幂法
    static Vector<Real> inverse_power(Matrix<Real> &A, Real lambda);   // 反幂法（半次迭代法）
    static int implicit_qr(Matrix<Real> &A, bool ifSymmetric = false); // 隐式 QR 方法
    static Matrix<Real> pass_Jacobi(Matrix<Real> &A, Real sigma);      // 过关 Jacobi 方法
    static Vector<Real> bisection(Matrix<Real> &A);                    // 二分法
};

// =====================================================

Matrix<Real> Solver::inverse(Matrix<Real> &A)
{
    int row = A.row();
    Matrix<Real> m(Matrix<Real>::unit(row));
    if (!A.square())
    {
        return m;
    }
    Matrix<Real> old(A);

    Real max;  // 最大值
    int index; // 最大值位置
    for (int i = 0; i < row; i++)
    {
        // find max
        max = 0;
        index = i;
        for (int j = i; j < row; j++)
        {
            Real tmp = fabs(old(j, i));
            if (tmp > max)
            {
                max = tmp;
                index = j;
            }
        }
        if (max == 0)
        {
            return m;
        }
        // change max
        for (int j = 0; j < row; j++)
        {
            Real tmp = old(i, j);
            old(i, j) = old(index, j);
            old(index, j) = tmp;

            tmp = m(i, j);
            m(i, j) = m(index, j);
            m(index, j) = tmp;
        }
        // start
        for (int j = 0; j < row; j++)
        {
            Real tmp = old(i, i);
            if (i == j)
            {
                for (int k = 0; k < row; k++)
                {
                    m(j, k) /= tmp;
                    old(j, k) /= tmp;
                }
                continue;
            }
            tmp = old(j, i) / old(i, i);
            for (int k = 0; k < row; k++)
            {
                m(j, k) -= tmp * m(i, k);
                old(j, k) -= tmp * old(i, k);
            }
        }
    }
    return m;
}

Real Solver::Householder(Vector<Real> &x, Vector<Real> &v)
{
    v = Vec::normalize(x, 0); // 防止溢出
    Real sumOfx = Vec::dot(v, v);
    Real sumOfx1 = sumOfx - v[0] * v[0];

    Real beta;
    if (sumOfx1 == 0)
    {
        beta = 0;
    }
    else
    {
        if (x[0] <= 0)
        {
            v[0] -= sqrt(sumOfx);
        }
        // 如果是正的分量，则用特殊算法减小舍入误差
        else
        {
            v[0] = -1 * sumOfx1 / (v[0] + sqrt(sumOfx));
        }
        // 对 beta 乘 v[0] * v[0] ，抵消规格化为 1 的影响
        beta = 2 * v[0] * v[0] / (sumOfx1 + v[0] * v[0]);
        v = v / v[0];
    }
    return beta;
}

Vector<Real> Solver::Givens(Vector<Real> &v, int i, int j)
{
    Real sqrt_x = sqrt(v[i] * v[i] + v[j] * v[j]);
    if (sqrt_x == 0)
    {
        return Vector<Real>::unit(2);
    }
    return Vector<Real>({v[i] / sqrt_x, v[j] / sqrt_x});
}

bool Solver::LU(Matrix<Real> &A)
{
    if (!A.square())
    {
        return false;
    }
    int n = A.row();
    for (int i = 0; i < n; i++)
    {
        Real key = A(i, i);
        // 无法进行高斯变换
        if (key == 0)
        {
            return false;
        }
        // 第 i 列高斯变换
        for (int j = i + 1; j < n; j++)
        {
            A(j, i) /= key;
            for (int k = i + 1; k < n; k++)
            {
                A(j, k) -= A(i, k) * A(j, i);
            }
        }
    }
    return true;
}

// 返回置换向量
void Solver::LU(Matrix<Real> &A, Perm &perm)
{
    if (!A.square())
    {
        return;
    }
    int n = A.row();
    Real max, tmp; // 记录最大值
    int index;     // 记录要交换的主元位置
    for (int i = 0; i < n; i++)
    {
        max = 0;
        index = i;
        // 选取列主元
        for (int p = i; p < n; p++)
        {
            tmp = fabs(A(p, i));
            if (tmp > max)
            {
                max = tmp;
                index = p;
            }
        }
        // 在这里奇异,说明 A 的秩为 i
        if (max == 0)
        {
            return;
        }
        // 交换主元
        for (int q = 0; q < n; q++)
        {
            tmp = A(i, q);
            A(i, q) = A(index, q);
            A(index, q) = tmp;
        }
        perm.exchange(i, index);
        // 正常的高斯消去法
        for (int j = i + 1; j < n; j++)
        {
            A(j, i) /= A(i, i);
            for (int k = i + 1; k < n; k++)
            {
                A(j, k) -= A(i, k) * A(j, i);
            }
        }
    }
}

// LL
bool Solver::Cholesky(Matrix<Real> &A)
{
    if (!A.square())
    {
        return false;
    }
    int n = A.row();
    Real sum = 0;
    for (int i = 0; i < n; i++)
    {
        sum = 0;
        for (int p = 0; p < i; p++)
        {
            sum += A(i, p) * A(i, p);
        }
        A(i, i) = sqrt(A(i, i) - sum);
        if (A(i, i) == 0)
        {
            return false;
        }
        for (int j = i + 1; j < n; j++)
        {
            sum = 0;
            for (int k = 0; k < i; k++)
            {
                sum += A(j, k) * A(i, k);
            }
            A(j, i) = (A(j, i) - sum) / A(i, i);
            A(i, j) = A(j, i);
        }
    }
    return true;
}

// LDL
bool Solver::Cholesky_D(Matrix<Real> &A)
{
    if (!A.square())
    {
        return false;
    }
    int n = A.row();
    Real sum = 0;
    for (int i = 0; i < n; i++)
    {
        sum = 0;
        for (int p = 0; p < i; p++)
        {
            sum += A(i, p) * A(i, p) * A(p, p);
        }
        A(i, i) -= sum;
        if (A(i, i) == 0)
        {
            return false;
        }
        for (int j = i + 1; j < n; j++)
        {
            sum = 0;
            for (int k = 0; k < i; k++)
            {
                sum += A(j, k) * A(i, k) * A(k, k);
            }
            A(j, i) -= sum;
            A(j, i) /= A(i, i);
            A(i, j) = A(j, i);
        }
    }
    return true;
}

void Solver::QR(Matrix<Real> &A, Vector<Real> &d)
{
    int column = A.col(), row = A.row();
    for (int i = 0; i < column; i++)
    {
        Vector<Real> x(row - i);
        Vector<Real> v(row - i);
        Real beta;
        for (int j = i; j < row; j++)
        {
            x[j - i] = A(j, i);
        }
        beta = Householder(x, v);
        Vector<Real> w(column - i);
        // get w = beta * AT * v
        for (int k = i; k < column; k++)
        {
            Real sum = 0;
            for (int p = i; p < row; p++)
            {
                sum += A(p, k) * v[p - i];
            }
            // note: it is w[k-i], not w[k]
            w[k - i] = beta * sum;
        }
        // get HA = A - v * wT
        for (int k = i; k < row; k++)
        {
            for (int p = i; p < column; p++)
            {
                if (p == i && k > i)
                {
                    A(k, p) = v[k - i];
                }
                else
                {
                    A(k, p) -= v[k - i] * w[p - i];
                }
            }
        }
        d[i] = beta;
    }
}

void Solver::rotate(Matrix<Real> &A, Matrix<Real> &Q, int p, int q)
{
    int n = A.row();
    Real cs[2];
    if (A(p, q) == 0)
    {
        cs[0] = 1;
        cs[1] = 0;
    }
    else
    {
        // t*t + 2*tau*t - 1 = 0
        Real tau = (A(q, q) - A(p, p)) / (2 * A(p, q));
        Real t = SIGN(tau);
        t /= (fabs(tau) + sqrt(1 + tau * tau));
        cs[0] = 1 / (sqrt(1 + t * t));
        cs[1] = t * cs[0];
    }
    // JT * A
    for (int j = 0; j < n; j++)
    {
        Real tmp;
        tmp = A(p, j);
        A(p, j) = tmp * cs[0] - A(q, j) * cs[1];
        A(q, j) = cs[1] * tmp + A(q, j) * cs[0];
    }
    // JT * A * J, note: here J is not a symmetical matrix
    // Q * J
    for (int j = 0; j < n; j++)
    {
        Real tmp;
        tmp = A(j, p);
        A(j, p) = tmp * cs[0] - A(j, q) * cs[1];
        A(j, q) = cs[1] * tmp + A(j, q) * cs[0];

        tmp = Q(j, p);
        Q(j, p) = tmp * cs[0] - Q(j, q) * cs[1];
        Q(j, q) = cs[1] * tmp + Q(j, q) * cs[0];
    }
}

void Solver::HAH(Matrix<Real> &A, int startOfRow, int endOfRow, int startOfCol, int endOfCol, Vector<Real> &v)
{
    // get houseHolder, beta, v
    // here householder will "not" change x, so we record the old length
    Real beta = v[0];
    v[0] = 1;
    Vector<Real> w(endOfCol - startOfCol);

    // caculate HA
    for (int j = startOfCol; j < endOfCol; j++)
    {
        Real sum = 0;
        // get w, (I - beta * v * vT) * A = A - v * wT
        // w = beta * AT * v
        for (int k = startOfRow; k < endOfRow; k++)
        {
            sum += A(k, j) * v[k - startOfRow];
        }
        w[j - startOfCol] = beta * sum;
    }
    // get A - v *wT
    for (int j = startOfRow; j < endOfRow; j++)
    {
        for (int k = startOfCol; k < endOfCol; k++)
        {
            A(j, k) -= v[j - startOfRow] * w[k - startOfCol];
        }
    }
    // caculate HAH, here we caculate H * (HA)T, so the form ti the same
    for (int j = startOfCol; j < endOfCol; j++)
    {
        Real sum = 0;
        for (int k = startOfRow; k < endOfRow; k++)
        {
            sum += A(j, k) * v[k - startOfRow];
        }
        w[j - startOfCol] = beta * sum;
    }
    // get A - v *wT
    for (int j = startOfRow; j < endOfRow; j++)
    {
        for (int k = startOfCol; k < endOfCol; k++)
        {
            A(k, j) -= v[j - startOfRow] * w[k - startOfCol];
        }
    }
}

void Solver::HAH(Matrix<Real> &A, int startOfRow, int endOfRow, Vector<Real> &v)
{
    // HAH = A - v *wT - w * vT
    Real beta = v[0];
    v[0] = 1;
    Vector<Real> u(endOfRow - startOfRow), w(endOfRow - startOfRow);

    // u = beta * A * v
    for (int i = startOfRow; i < endOfRow; i++)
    {
        Real sum = 0;
        for (int j = startOfRow; j < endOfRow; j++)
        {
            sum += A(i, j) * v[j - startOfRow];
        }
        u[i - startOfRow] = beta * sum;
    }
    // w = u - 1/2 * beta * (vT *u) * v
    w = u - v * (v * u) * beta * 0.5;
    for (int i = startOfRow; i < endOfRow; i++)
    {
        for (int j = startOfRow; j < endOfRow; j++)
        {
            A(i, j) -= v[i - startOfRow] * w[j - startOfRow] + v[j - startOfRow] * w[i - startOfRow];
        }
    }
}

void Solver::Hessenberg(Matrix<Real> &A, Matrix<Real> &Q, bool ifSymmetric)
{
    int n = A.row();
    for (int i = 0; i < n - 2; i++)
    {
        Vector<Real> x = A(i).range(i + 1, n), v(n - i - 1);
        Real beta = Householder(x, v);
        v[0] = beta;
        ifSymmetric ? HAH(A, i + 1, n, v) : HAH(A, i + 1, n, i, n, v);
        Q.insert(i + 1, i, v);
    }
}

void Solver::Real_QR(Matrix<Real> &A, int l, int m)
{
    int n = A.row();
    int _m = n - 1;
    Real s = A(_m - 1, _m - 1) + A(_m, _m);
    Real t = A(_m - 1, _m - 1) * A(_m, _m) - A(_m - 1, _m) * A(_m, _m - 1);
    Vector<Real> x(3);
    x[0] = A(0, 0) * A(0, 0) + A(0, 1) * A(1, 0) - s * A(0, 0) + t;
    x[1] = A(1, 0) * (A(0, 0) + A(1, 1) - s);
    x[2] = A(1, 0) * A(2, 1);
    // here we ignore the position that n<3!!!!!!!!!!!!!
    // get P0 * H * P0

    Vector<Real> v(3);
    Real beta;
    beta = Householder(x, v);
    v[0] = beta;
    HAH(A, l, l + 3, l, n, v);

    x[0] = A(l + 1, l);
    x[1] = A(l + 2, l);
    if (l + 3 < n - m)
    {
        x[2] = A(l + 3, l);
        for (int i = l + 1; i < n - m - 2; i++)
        {
            // householder
            beta = Householder(x, v);
            v[0] = beta;
            HAH(A, i, i + 3, i - 1, n, v);

            x[0] = A(i + 1, i);
            x[1] = A(i + 2, i);
            if (i < n - m - 3)
            {
                x[2] = A(i + 3, i);
            }
        }
    }
    Vector<Real> y({x[0], x[1]});
    Vector<Real> v2(2);
    beta = Householder(y, v2);
    v2[0] = beta;
    HAH(A, n - m - 2, n - m, n - m - 3, n, v2);
}

void Solver::Wilkinson_QR(Matrix<Real> &A, int l, int m)
{
    int n = A.row();
    int last = n - m - 1;
    Real d = (A(last - 1, last - 1) - A(last, last)) / 2;
    int sign = SIGN(d);
    Real u = A(last, last) - A(last, last - 1) * A(last, last - 1) / (d + sign * sqrt(d * d + A(last, last - 1) * A(last, last - 1)));
    Vector<Real> x({A(l, l) - u, A(l + 1, l)});
    Vector<Real> cs(2);
    for (int i = l; i < last; i++)
    {
        cs = Givens(x);
        // GT
        for (int j = 0; j < n; j++)
        {
            Real tmp = 0;
            tmp = A(i, j);
            A(i, j) = tmp * cs[0] + A(i + 1, j) * cs[1];
            A(i + 1, j) = -1 * cs[1] * tmp + A(i + 1, j) * cs[0];
        }
        // GTG^T, note: here G is not a symmetical matrix
        for (int j = 0; j < n; j++)
        {
            Real tmp = 0;
            tmp = A(j, i);
            A(j, i) = tmp * cs[0] + A(j, i + 1) * cs[1];
            A(j, i + 1) = -1 * cs[1] * tmp + A(j, i + 1) * cs[0];
        }
        if (i < last - 1)
        {
            x[0] = A(i + 1, i);
            x[1] = A(i + 2, i);
        }
    }
}

// ===================================================================================================

Vector<Real> Solver::solve_upper(Matrix<Real> &A, Vector<Real> &b, bool unit)
{
    int n = b.size();
    Vector<Real> x(n);
    Real sum;
    for (int i = n - 1; i >= 0; i--)
    {
        sum = 0;
        for (int j = i + 1; j < n; j++)
        {
            sum += A(i, j) * x[j];
        }
        x[i] = (b[i] - sum) / (unit ? 1 : A(i, i));
    }
    return x;
}

Vector<Real> Solver::solve_lower(Matrix<Real> &A, Vector<Real> &b, bool unit)
{
    int n = b.size();
    Vector<Real> x(n);
    Real sum;
    x[0] = b[0] / (unit ? 1 : A(0, 0));
    for (int i = 1; i < n; i++)
    {
        sum = 0;
        for (int j = 0; j < i; j++)
        {
            sum += A(i, j) * x[j];
        }
        x[i] = (b[i] - sum) / (unit ? 1 : A(i, i));
    }
    return x;
}

Vector<Real> Solver::solve(Matrix<Real> &A, Vector<Real> &b, bool choose)
{
    int n = b.size();
    Vector<Real> x(n);
    if (choose)
    {
        Perm perm(n);
        LU(A, perm);
        b = perm * b;
    }
    else
    {
        LU(A);
    }
    // 单位下三角阵
    x = solve_lower(A, b, true);
    x = solve_upper(A, x);
    return x;
}

Vector<Real> Solver::cholesky(Matrix<Real> &A, Vector<Real> &b, bool improve)
{
    int n = b.size();
    Vector<Real> x(n);
    if (improve)
    {
        Cholesky_D(A);
        // 将对角元乘到上三角部分
        for (int i = 0; i < n; i++)
        {
            for (int j = i + 1; j < n; j++)
            {
                A(i, j) *= A(i, i);
            }
        }
        x = solve_lower(A, b, true);
        x = solve_upper(A, x);
    }
    else
    {
        Cholesky(A);
        x = solve_lower(A, b);
        x = solve_upper(A, x);
    }
    return x;
}

Vector<Real> Solver::qr(Matrix<Real> &A, Vector<Real> &b)
{
    int row = A.row(), col = A.col();
    Vector<Real> d(col);
    QR(A, d);

    // Householder 变换作用于 b
    // Hb = b - beta * v * v^T * b = b - (beta * v^T * b) * v
    for (int i = 0; i < col; i++)
    {
        Real sum = b[i];
        for (int j = i + 1; j < row; j++)
        {
            sum += A(j, i) * b[j];
        }
        sum *= d[i], b[i] -= sum;
        for (int j = i + 1; j < row; j++)
        {
            b[j] -= A(j, i) * sum;
        }
    }
    // 求解上半部分的线性方程组
    Vector<Real> b1(b.range(0, col));
    Vector<Real> x(solve_upper(A, b1));
    return x;
}

Vector<Real> Solver::Jacobi(Matrix<Real> &A, Vector<Real> &b)
{
    int n = A.order();
    // x = D1(L+U) * x + D1 * g
    Vector<Real> x(n, 1), y(n); // let x be a all one vector
    Real sum, dx;
    do
    {
        dx = 0;
        // 分步迭代过程
        for (int i = 0; i < n; i++)
        {
            sum = 0;
            for (int j = 0; j < n; j++)
            {
                if (i != j)
                {
                    sum -= A(i, j) * x[j];
                }
            }
            y[i] = (sum + b[i]) / A(i, i);
        }
        dx = Vec::norm(y - x, 0);
        x = y;
    } while (fabs(dx) >= ACCURACY);
    return x;
}

Vector<Real> Solver::G_S(Matrix<Real> &A, Vector<Real> &b)
{
    int n = A.order();
    Real sum, dx;
    Vector<Real> x(n, 1);

    do
    {
        dx = 0;
        // 分部迭代过程
        for (int i = 0; i < n; i++)
        {
            sum = 0;
            for (int j = 0; j < n; j++)
            {
                if (j != i)
                {
                    sum -= A(i, j) * x[j];
                }
            }
            Real tmp = x[i];
            x[i] = (sum + b[i]) / A(i, i);
            dx = MAX(fabs(tmp - x[i]), dx);
        }
    } while (dx >= ACCURACY);
    return x;
}

Vector<Real> Solver::SOR(Matrix<Real> &A, Vector<Real> &b, Real w)
{
    // x = Lw * x + w(D-wL)1 * b
    int n = A.order();
    Vector<Real> x(n, 1);
    Real sum, dx;

    do
    {
        dx = 0;
        // 分部迭代过程
        for (int i = 0; i < n; i++)
        {
            sum = 0;
            for (int j = 0; j < n; j++)
            {
                if (j != i)
                {
                    sum -= A(i, j) * x[j];
                }
            }
            Real tmp = x[i];
            x[i] = (1 - w) * x[i] + w * (sum + b[i]) / A(i, i);
            dx = MAX(fabs(tmp - x[i]), dx);
        }
    } while (dx >= ACCURACY);
    return x;
}

Real Solver::kappa(Matrix<Real> &A)
{
    int n = A.order();
    Real A_K = Mat::norm(A, 0); // A's condition number
    Real A1_K = 0;              // A-1's condition number
    Matrix<Real> AT(Mat::transpose(A));
    Matrix<Real> _A(A);
    Vector<Real> x(n, 1.0 / n), w(n), v(n), z(n);

    // 分别对 AT 和 A 做分解
    LU(AT);
    LU(A);

    // 判断是否终止
    bool yes = false;
    int j = 0;

    while (!yes)
    {
        // get w
        w = solve_lower(AT, x);
        w = solve_upper(AT, w);
        v = static_cast<Vector<Real>>(Vec::sgn(w));
        // get z
        z = solve_lower(A, v);
        z = solve_upper(A, z);
        // compare infinite z and zT * x
        if (Vec::norm(z, 0) <= Vec::dot(z, x))
            break;
        x = Vector<Real>::unit(n, Vec::max_index_of(z));
    }
    // A1_K = one w
    A1_K = Vec::norm(w, 1);
    return A_K * A1_K;
}

Real Solver::caculate_analysis(Matrix<Real> &A, Vector<Real> &x, Vector<Real> &b)
{
    int n = A.order();
    Vector<Real> r(A * x - b);
    Real kr, kb, K;
    kr = Vec::norm(r, 0);
    kb = Vec::norm(b, 0);
    K = kappa(A);
    return K * kr / kb;
}

Vector<Real> Solver::cg(Matrix<Real> &A, Vector<Real> &b)
{
    int row = A.row(), column = A.col(), n = ITERATION;
    Vector<Real> x(row, 1), p(row), r(row);
    p = b - A * x, r = p;                     // get p_0 = r_0
    Real rr = Vec::dot(r, r), nrr = rr, beta; // get r0 * r0, next r0 * r0
    Real alpha = rr / Vec::dot(p, A * p);     // get alpha_0 = r0T * r0 / p0T * A * p0
    x = p * alpha + x;                        // get x_0

    while (sqrt(rr) > ACCURACY && n > 0)
    {
        r = r - A * p * alpha;           // get r = rk-1 - alphak-1 * A * pk-1
        nrr = Vec::dot(r, r);            // get next r * r, nrr = r * r
        beta = nrr / rr;                 // get beta = rkT * rk / rk-1T * rk-1
        rr = nrr;                        // get p = rk + betak-1 * pk-1
        p = r + p * beta;                // get alpha
        alpha = rr / Vec::dot(p, A * p); // get x
        x = p * alpha + x;
        n--;
    }
    return x;
}

Real Solver::power(Matrix<Real> &A, Real mu)
{
    int n = A.order();
    Vector<Real> x(n, 1);
    Matrix<Real> _A(n, n);
    _A = A - x * mu;

    Real dp = 1, p, q = 0;
    Vector<Real> y(n);
    while (dp > ACCURACY)
    {
        y = _A * x;
        p = y[Vec::max_index_of(y)];
        x = y / p;
        dp = fabs(p - q);
        q = p;
    }
    return p + mu;
}

Vector<Real> Solver::inverse_power(Matrix<Real> &A, Real lambda)
{
    int n = A.order();
    Vector<Real> e(n, 1);
    Matrix<Real> _A(n, n);
    _A = A - e * lambda;
    // z_0 = P^T * L * e, U * v_1 = e
    // z = PermOp::inverse(p) * (lower(_A) + z) * z;

    Perm p(n);
    LU(_A, p);
    e = solve_upper(_A, e);
    return Vec::normalize(e, 0);
}

int Solver::implicit_qr(Matrix<Real> &A, bool ifSymmetric)
{
    int l, m, n = A.order(), count = 0;
    Matrix<Real> Q(n, n);
    Hessenberg(A, Q, ifSymmetric);
    while (count < ITERATION)
    {
        // Zero
        for (int i = 1; i < n; i++)
        {
            if (fabs(A(i, i - 1)) <= ((fabs(A(i - 1, i - 1)) + fabs(A(i, i))) * ACCURACY))
            {
                A(i, i - 1) = 0;
                if (ifSymmetric)
                {
                    A(i - 1, i) = 0;
                }
            }
        }
        for (int i = 2; i < n; i++)
        {
            for (int j = 0; j < i - 1; j++)
            {
                A(i, j) = 0;
                if (ifSymmetric)
                {
                    A(j, i) = 0;
                }
            }
        }
        // divide to three part
        m = 0;
        Real last = A(n - 1, n - 2);
        for (int i = n - 2; i > 0; i--)
        {
            if (A(i, i - 1) == 0 || last == 0)
            {
                m++;
            }
            else
            {
                break;
            }
            last = A(i, i - 1);
            if (i - 1 == 0)
            {
                m = n;
            }
        }
        if (m == n)
        {
            break;
        }
        else
        {
            l = 0;
            for (int i = n - m - 1; i > 0; i--)
            {
                if (A(i, i - 1) != 0)
                {
                    l++;
                }
            }
            l++;
            l = n - m - l;
        }
        // Real QR
        if (ifSymmetric)
        {
            Wilkinson_QR(A, l, m);
        }
        else
        {
            Real_QR(A, l, m);
        }
        count++;
    }
    return count;
}

Matrix<Real> Solver::pass_Jacobi(Matrix<Real> &A, Real sigma)
{
    int n = A.order();
    Matrix<Real> Q(Matrix<Real>::unit(n));
    if (sigma < n)
    {
        return Q;
    }
    Real delta = Mat::norm(A, 'e');
    bool check;
    int count = 0;
    while (delta > ACCURACY && count < ITERATION)
    {
        check = true;
        for (int i = 0; i < n; i++)
        {
            for (int j = i + 1; j < n; j++)
            {
                if (fabs(A(i, j)) > delta)
                {
                    rotate(A, Q, i, j);
                    check = false;
                }
            }
        }
        if (check)
        {
            delta /= sigma;
        }
        count++;
    }
    return Q;
}

Vector<Real> Solver::bisection(Matrix<Real> &A)
{
    int n = A.order();
    Vector<Real> lambda(n), alpha(n), beta(n);
    Real u = Mat::norm(A, 0);
    Real l = -1 * u;
    for (int i = 0; i < n; i++)
    {
        alpha[i] = A(i, i);
        if (i > 0)
        {
            beta[i] = A(i, i - 1);
        }
    }
    for (int i = 0; i < n; i++)
    {
        Real _u = u, _l = l;
        Real nowE = u - l, r = (u + l) / 2;
        while (nowE > ACCURACY)
        {
            // compare i+1 and sign reverse number
            int s = 0;
            Real q = alpha[0] - r;
            for (int j = 0; j < n; j++)
            {
                if (q < 0)
                {
                    s++;
                }
                if (j < n - 1)
                {
                    if (q == 0)
                    {
                        // replace q with a little positive number
                        q = fabs(beta[j + 1]) * ACCURACY;
                    }
                    q = alpha[j + 1] - r - beta[j + 1] * beta[j + 1] / q;
                }
            }
            if (s >= (i + 1))
            {
                _u = r;
            }
            else
            {
                _l = r;
            }
            r = (_u + _l) / 2;
            nowE = _u - _l;
        }
        lambda[i] = r;
    }
    return lambda;
}